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Egypt has the highest prevalence of hepatitis C virus (HCV) infection worldwide with a frequency of 15%. More than 
90% of these infections are due to genotype 4, and the subtype 4a (HCV-4a) predominates. Moreover, due to the 
increased mobility of people, HCV-4a has recently spread to several European countries. The protease domain of the 
HCV nonstructural protein 3 (NS3) has been targeted for inhibition by several drugs. This approach has had marked 
success in inhibiting genotype 1 (HCV-1), the predominant genotype in the USA, Europe, and Japan. However, HCV-4a 
was found to resist inhibition by a number of these drugs, and little progress has been made to understand the structural 
basis of its drug resistivity. As a step forward, we sequenced the NS3 HCV-4a protease gene (strain ED43) and subse- 
quently built a 3D structural model threaded through a template crystal structure of HCV-lb NS3 protease. The model 
protease, HCV-4a, shares 83% sequence identity with the template protease, HCV-lb, and has nearly identical rigid 
structural features. Molecular dynamics simulations predict similar overall dynamics of the two proteases. However, local 
dynamics and 4D analysis of the interactions between the catalytic triad residues (His57, Asp81, and Serl39) indicate 
conformational instability of the catalytic site in HCV-4a NS3 protease. These results suggest that the divergent dynam- 
ics behavior, more than the rigid structure, could be related to the altered catalytic activity and drug resistivity seen in 
HCV-4a. 
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Introduction 

Hepatitis C virus (HCV) is a global health concern. 
Chronic infection of HCV is a common and leading 
cause for both cirrhosis and hepatocellular carcinoma 
(Andrade et al., 2009). Approximately 3% of the world's 
population, or roughly 170 million people, are currently 
affected by this disease (Levanchy, 2009). Compara- 
tively, the frequency within Egypt is much higher with 
15%, or nearly 13 million people, testing HCV seroposi- 
tive (Miller & Abu-Raddad, 2010). Roughly 90% of the 
aforementioned Egyptians are carriers of HCV genotype 
4, and the subtype 4a (HCV-4a) predominates (Khattab 
et al., 2011; Nguyen & Keeffe, 2005). 

Genotype 1 is the common variant of HCV through- 
out the USA, Europe, and Japan, and has thus become 
the focus of much interest and research (Ali, Ahmed, & 
Idrees, 2010). While no vaccine is available and current 
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therapies had met with limited success, the management 
and treatment of infections arising from genotype 1 has 
advanced considerably (Chatel-Chaix, Baril, & Lamarre 
2010; Kwo & Vinayek, 2011). In contrast, however, 
genotype 4 has not undergone adequate scrutiny and as a 
result, the targeted drug development has stagnated 
(Kamal & Nasser, 2008). 

Due to its importance in the replication cycle of 
HCV, the serine protease domain of nonstructural protein 
3 (NS3) has been an attractive target for the development 
of effective inhibitors (Heintges, Enche, Putlitz, & 
Wands, 2001). The NS3 protease cleaves four down- 
stream sites in the polyprotein and is characterized as a 
serine protease with a chymotrypsin-like fold, which is 
activated by the NS4A cofactor (Du, Hou, Guan, Tong, 
& Wang, 2002). Similar to chymotrypsin, the catalytic 
triad of HCV NS3 protease is comprised of the three 
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essential residues, histidine, aspartic acid, and serine, 
numbered from the N-terminus of NS3, 57, 81, and 139, 
respectively (Lin, 2006). 

Together, these three residues carry out general 
acid-base catalysis on target peptides. Throughout the 
catalytic mechanism, two tetrahedral intermediates are 
formed. The serine performs a nucleophilic attack on a 
carbonyl of the substrate and the histidine fulfills the 
crucial roles as both acid and base that allow the 
catalysis to progress (Hedstrom, 2002). The aspartic acid 
stabilizes histidine via hydrogen bonding that also raises 
the pKa value of the histidine which is essential in 
catalysis (Fersht & Sperling, 1973; Hedstrom, 2002). 

Other interactions have also been shown to affect 
catalysis. Zinc plays an important role in the structural 
stability of NS3 protease by enthalpically disfavoring 
protein denaturation (Abian, Vega, Neira, & Valazquez- 
Campoy, 2010). Additionally, a bound peptide cofactor 
(NS4A) increases the protease activity by nearly 
1000-fold (Sardana, Blue, Zugay-Murphy, Sardana, & 
Kuo, 1999). 

A several-fold decrease in the catalytic efficiency of 
HCV-4a NS3 protease has been reported relative to that 
of HCV-lb (Franco, Clotet, & Martinez, 2008). Several 
NS3 protease inhibitors, which were designed to interfere 
with the catalytic triad (e.g. Telaprevir and Boceprevir) 
have shown promising results with genotype 1, but not 
with genotype 4 (Chatel-Chaix et al., 2010; Njoroge, 
Chen, Shih, & Piwinski 2008; White et al, 2010). Very 
little has been reported on the basis of HCV-4's 
drug resistance, mainly due to the lack of available 3D 
structural information. 

Molecular dynamics simulations have been useful in 
the study of the NS3 protease domain of HCV. These 
simulations elucidated the interaction between the NS4A 
cofactor and the NS3 protease of genotype lb (Zhu & 
Briggs, 2011). They have also granted an insight into the 
effect of R155K, A156V, and D168A mutations on the 
resistance of HCV-lb to the protease inhibitors, ITMN- 
191 and TMC435 (Pan, Xue, Zhang, Liu, & Yao, 2012; 
Xue, Pan, Yang, Liu, & Yao, 2012). 

In the current study, we resequenced the gene coding 
for the HCV-4a (strain ED43) NS3 protease domain and 
subsequently built a novel 3D structural model including 
a bound NS4A cofactor and a zinc ion. The model was 
threaded through the X-ray crystal structure of HCV-lb 
protease (PDB: ldy8, Di Marco et al., 2000). We per- 
formed molecular dynamics simulations, on both the 
threading model and the template structure, to investigate 
the effect of the sequence variability on their respective 
dynamic behaviors. 

Our modeling and molecular dynamics simulations 
results show that both proteases share very similar rigid 
and overall dynamic features. Conversely, both proteases 
exhibit significantly different local dynamics and distance 



distribution profiles, both in peak values and broadness, 
at the catalytic triad. Our data suggest that genotype- 
based structural dynamics could play a significant role in 
the stability of the catalytic triad, and thereby, in drug 
response among the HCV genotypes. 

Materials and methods 

DNA sequencing 

Isolate ED43 of the HCV-4a genome, inserted in PUC19 
vector, was kindly provided by the lab of Dr. Richard 
M. Elliott (University of Glasgow, Glasgow, Scotland, 
UK). The complete nucleotide sequence of HCV-ED43, 
which is 9355 nucleotides long, is deposited in the 
EMBL database under accession number Y11604 
(Chamberlain, Adams, Saeed, Simmonds, & Elliott, 
1997). The DNA sequence coding for the NS3 proteins, 
including the protease, was cloned into pQE-30 vector 
using pQE-30 primers type III/IV as previously described 
(Sabet, Al-Sherbiny, Ibrahim, & Hagan, 2009). 

To confirm the sequence results, both the template 
and cloned DNA were sequenced. The two sequencings 
were done at MBSU DNA sequencing services at the 
Institute of Biomedical and Life Sciences, University of 
Glasgow, Scotland, UK. These sequencings were per- 
formed using automated sequencer, ABI PRISM model 
377 version 3.3.1. The results obtained from the two 
sequences were aligned to the original ED43 published 
sequence Y11604. Three sites T54, 1134, and R150 in 
the published sequence were found to be 154, T134, and 
A150, respectively. The gene coding for the NS4A 
cofactor was sequenced separately using ABI PRISM 
model 3730XL analyzer and the result was identical to 
the published sequence Y11604. 

3D structure prediction and validation 

The 3D structure of HCV-4a NS3 protease was predicted 
by threading its amino acid sequence through the X-ray 
crystal structure of HCV-lb NS3 protease(ldy8) via the 
threading program LOOPP (Di Marco et al, 2000; Mel- 
ler & Elber, 2001). LOOPP is a fold recognition program 
that generates atomic coordinates of a sample molecule 
based on an alignment with a homologous template 
structure. By integrating the results from direct sequence 
alignment, sequence profile, threading, secondary struc- 
ture, and exposed surface area prediction, the LOOPP 
builds main-chain and all-atom models. A nearly identi- 
cal model (RMSD 0.1 56 A) was also obtained via 
homology modeling using the SWISS-MODEL Work- 
space (Arnold, Bordoli, Kopp, & Schwede, 2006; Guex, 
& Peitsch, 1997; Kiefer, Arnold, Kunzli, Bordoli, & 
Schwede, 2009; Peitsch, 1995; Schwede, Kopp, Guex, & 
Peitsch, 2003). 

To build the NS4A cofactor, we superposed the model 
structure onto the template structure (ldy8, RMSD 0.3 A) 
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and built the sequence of the NS4A cofactor for the 
model based on the corresponding coordinates found in 
the template crystal structure. Similarly, a single zinc ion 
was manually docked at the cysteine triad C97, C99, and 
CI 45 into the model guided by the corresponding posi- 
tion in another structure of the template protein (ldxp) in 
which zinc is present (Di Marco et al., 2000). With the 
cofactor and zinc bound, the model was energy mini- 
mized using the CCP4 program suite (Potterton, Briggs, 
Turkenburg, & Dodson, 2003; Winn et el, 2011) and the 
GROMOS96 program, an implementation of the Swiss- 
pdb viewer (Van Gunsteren et al., 1996). 

The final model was validated using the NIH MB1 
Laboratory for Structural Genomics and Proteomics 
Structural Analysis and Verification Server. This server 
utilizes five programs (Procheck, What Check, ERRAT, 
Verify_3D, and Prove) to analyze the stereochemical 
parameters and the quality of the model (Colovos & 
Yeates, 1993; Hooft, Vriend, Sander, & Abola, 1996; 
Laskowski, MacArthur, Moss, & Thorton, 1993; Luthy, 
Bowie, & Eisenberg, 1992; Pontius, Richelle, & Wodak, 
1996). 

Additionally, CCP4 programs suite 6.0 was used for 
the calculation of a Ramachandran plot (Ramachandran, 
Ramakrishnan, & Sasisekharan, 1963), structure superpo- 
sition, and rmsd value calculation in addition to the eval- 
uation of the stereochemistry (Winn et al, 2011). 

Molecular dynamics simulation 

The molecular dynamics simulation (MD) was performed 
using NAMD 2.9 under the CHARMM27 force field for 
proteins (MacKerell et al., 1998; MacKerell, Feig, & 
Brooks, 2004; Phillips et al, 2005). Initially, the 3D 
structure was solvated using the solvation tool in VMD 
(Humphrey, Dalke, & Schulten, 1996). The TIP3P model 
was used for the water molecules (Jorgensen, Chandrase- 
khar, Madura, Impey, & Klein, 1983). Lengevin dynam- 
ics for all nonhydrogen atoms with a damping coefficient 
of 1 ps -1 was used in maintaining a constant temperature 
of 310K throughout the system. A constant pressure of 
1 atm was maintained using a Nose-Hoover Langevin 
piston with a period of 1 00 fs and damping timescale of 
50 fs (Feller, Zhang, Pastor, & Brooks, 1995). 

Periodic boundary conditions were used on a 61 A 
cubic box with the long-range electrostatics calculated 
using the particle-mesh Ewald method with a grid point 
density of 0.92 A -1 . This process ensured that adjacent 
copies of the protease were never close enough for 
short-range interaction. A cut-off of 10 A for van der 
Waals interactions and a switching distance of 8 A were 
found to give convergent results, thus used for production 
runs. The solvation box was neutralized, using VMD's 
Autoionize plugin version 1.3, with sodium chloride 
placed at distances greater than 5 A from the protease. 



A time step of 1 fs was used in order to resolve the 
hydrogen motion of water. The initial structure was first 
subjected to two rounds of an 800 cycle conjugate gradi- 
ent energy minimization followed by 50 ps of MD simu- 
lation at 100K. The system was then heated up in 
increments of 25 K with 50 ps of MD simulation at each 
temperature increment until the desired temperature of 
310 K was established. The system was then simulated 
for 30 ns. Time frames used for the measurements within 
the protease were only done for frames where the prote- 
ase had equilibrated: 6-25 ns. The equilibrium state of 
the protease was determined by the RMSD of the entire 
protein's backbone. 

Multiple copies of each protease, which included the 
cofactor and a zinc ion (nonbonded), were run with dif- 
ferent initial conditions to ensure that the results were 
well converged. All RMSD data presented are averaged 
over for at least three distinct runs. Distance distributions 
are for a single run that was typical of the set of runs to 
ensure the clarity of interpretation. 

Results and discussion 

The NS3 protease domain of HCV is regarded as a 
chymotrypsin-like serine protease (Tong, 2002) whose 
catalytic mechanism is well established (Hedstrom, 
2002). In summary, a charge relay system is formed in 
which the carboxylic group of D81 forms a hydrogen 
bond with N<51 of H57 (Hedstrom, 2002). This event 
increases the pKa of the histidine side chain from 7 to 
about 12 (Fersht & Sperling, 1973; Lin, Cassidy, & Frey, 
1998). Consequently, H57 deprotonates the hydroxyl 
group of the SI 39 side chain and a proton shuttles to 
Ne2 of His57 (Hedstrom, 2002). The Oy of S139 then 
nucleophilically attacks the carbonyl carbon of a sub- 
strate's scissile bond resulting in the formation of an 
oxyanion-containing tetrahedral intermediate (Hedstrom, 
2002; Hirohara, Bender, & Stark, 1974; Raney, Sharma, 
Moustafa, & Cameron, 2010; Tong, 2002). At this point, 
the protonated H57 acts as a general acid assisting in the 
collapse of the tetrahedral intermediate and the cleavage 
of the substrate (Hedstrom, 2002; Raney et al., 2010). 
Given the function of the NS3 protease within the viral 
replication process, the catalytic triad has become an 
attractive target for rational drug design. 

The HCV-4a NS3 protease structure model super- 
poses very well on the threading template structure 
(ldy8), and the two share 83% sequence identity 
(Figure 1(d)) along with nearly identical structural fea- 
tures (Figure 1(a)). RMSD in back-bone positions 
between the two proteases is about 0.3 A. None of the 
174-threaded amino acids falls within the disallowed 
Ramachandran area, and no steric clashes or stereochem- 
ical outliers was detected (see Materials and Methods). 
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Figure 1. Comparison between the threading model of HCV-4a and the crystal structure of HCV-lb NS3 proteases, (a) The model 
structure of HCV-4a, shown in blue, is superimposed onto the template structure of HCV-lb (PDB: ldy8), shown in green. Sequence 
variability sites are shown in red. The tan box highlights the catalytic triad (H57, D81 and S139). A nearby mutation T54I is shown 
in red sticks, (b) Same as (a) with the addition of the cofactor (NS4A), in pink, and an example inhibitor shown for orientation 
purposes (OF7, ldy8), in orange, (c) Residue-averaged RMSD of Ca atoms for HCV-4a protease model (blue) and HCV-lb protease 
crystal structure (green) during the course of simulation, (d) The amino acid sequences of HCV-4a and HCV-lb NS3 proteases, as 
well as the cofactor are shown with the same colors used in (a) and (b). Sequence variability is indicated with a red "X". The 
position of T54I is highlighted by a bold underlined red "X". 



The three catalytic residues H57, D81, and S139 are 
located in a crevice between the two protease (3-barrels 
as shown in Figure 1(a) (Barbato et al., 1999; McCoy, 
Senior, Gesell, Ramanathan, & Wyss, 2001; Yan et al., 
1998). The active site is nonpolar and shallow (Hed- 
strom, 2002). The central region of NS4A is buried 
almost completely inside the NS3 protease and serves as 
a cofactor for proper folding of the protease (Figure 1 
(b)) (Barbato et al, 1999). The rigid structures indicate 
that access to the active site is nearly identical in both 
the model and template (Figure 1(a) and (b)). In addi- 
tion, molecular dynamics simulations predict that both 
HCV-4a and HCV-lb proteases share similar average 
RMSD in the Ca positions, ~1.3A at equilibrium (Fig- 
ure 1(c)). Locally, molecular dynamics simulations 
revealed divergent dynamics behavior and distance distri- 
bution profiles within the catalytic triad region between 
HCV-4a and HCV-lb proteases (Figures 2 and 3). These 
dynamic differences are likely to contribute to the altered 
activity and drug resistivity seen in HCV-4a. 



The dynamical behavior within the catalytic triad 

The alpha carbons (Ca) of the catalytic residues exhi- 
bit similar dynamics throughout the simulations with 
those of the HCV-4a model showing slightly higher 
RMSD, especially for H57 and D81 (Figure 2(d-f)). 
The SI 39 side chain also presents a very similar 
dynamic behavior in both the model (HCV-4a) and 
template (HCV-lb) (Figure 2(c)). These data are 
expected since the behavior of Ca atoms, as part of 
the peptide backbone, would conform to the dynamics 
seen in Figure 1(c). In addition, the serine presents a 
relatively small side chain with a movement sterically 
encumbered by the nearby residue. 

However, in the HCV-4a model, H57 and D81 (as 
entire residues) demonstrate dynamics behavior divergent 
from that of the template HCV-lb (Figure 2(a) and (b)). 
The RMSD of H57 in the HCV-4a model varies from 
that predicted in the HCV-lb template by up to 1.5 A at 
certain points during the simulation (Figure 2(a)). Like- 
wise, the RMSD of D81 in the HCV-4a model diverges 
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Figure 2. Comparison of the dynamics behavior of the catalytic triad residues between the threading model of HCV-4a (blue) and 
the crystal structure of HCV-lb (green) proteases. RMSD values for each catalytic residue are shown for the entire residue (a, b, c) 
and the corresponding alpha carbons (d, e, f). 



from that of HCV-lb template by nearly 1 A (Figure 2 
(b)). It is possible that the instability of D81, itself 
fulfilling a stabilizing role for H57, may have deleterious 
consequences upon the ability of H57 to adequately 
function as an acid-base. Given the crucial role of H57 
in the proteolytic mechanism, our data provide a likely 
explanation for the reported reduced efficiency of the 



HCV-4a NS3 protease compared with that of HCV-lb 
(Franco et al, 2008). 

4D simulation of the interactions between the catalytic 
residues 

In order to assess the potential impact of the aforemen- 
tioned differences, we investigated the positional dynam- 
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Figure 3. Dynamics behavior within the catalytic triad site of the threading model (HCV-4a, blue) and the template (HCV-lb, 
green) proteases. The distance distribution profiles between 081 of residues D81 and N51 of H57 (a) and between Oy of residue 
S139 and Ne2 of residue H57 (b), during the course of simulation. Orange arrows indicate the selected distances in the rigid 
structures of both the model (c) and the template (d). Residues shown in red indicate sequence variability. 



ics of the catalytic triad residues during the course of 
simulation and used the distance distribution profiles of 
catalytically relevant distances as indicators of the 4D 
variations. 

The distance distribution profiles between Ne2 of 
H57 and Oy of S139, as well as between N51 of H57 
and 081 of D81, of the template (HCV-lb) and model 
(HCV-4a) vary widely in both peak value and breadth. 

In the template structure (HCV-lb), the distance 
between N81 of H57 and 051 of D81 exhibits a sharp 
distribution with a peak value around 3 A. In the model 
(HCV-4a), however, the corresponding distance distribu- 
tion is bimodal, much broader, and distributed around 
^7 and 10 A (Figure 3(a)). Similar distributions were 
obtained between N81 of H57 and 052 of D81 (see Sup- 
plementary data). Furthermore, the distance distribution 
between Oy of SI 39 and Ne2 of H57 in the template 
(HCV-lb) shows a peak value at around 4 A, while the 
corresponding distribution in the model (HCV-4a) is 
broader and exhibits a shifted peak value at around 8 A 
(Figure 3(b)). 



Together, these data indicate that in the model (HCV- 
4a), H57 spends less time within a probable hydrogen 
bonding distance to both SI 39 and D81. Thus, H57 is 
less likely to act as an efficient general acid-base, as 
explained previously. This is consistent with the observa- 
tion that the catalytic activity of HCV-4a NS3 protease is 
several folds of magnitude less than that of HCV-lb 
(Franco et al., 2008). 

It is important to note that the predicted divergent 
dynamics behavior (Figure 3(a) and (b)) is totally hidden 
by the apparent similarity seen in the catalytic site of the 
rigid structures (Figure 3(c) and (d)). This highlights the 
importance of utilizing molecular dynamics as a method 
of future investigations into protease activity. 

T54I mutation may affect the dynamics of His57 

HCV NS3 protease inhibitors such as Telaprevir and 
Boceprevir are generally designed to interfere with the 
catalytic triad (Kwo & Vinayek, 2011; Raney et al., 
2010). Precedence for drug resistance due to variability 
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Figure 4. Dynamics behavior near the catalytic H57 in the model (HCV-4a, blue) and the template (HCV-lb, green) proteases (b,c). 
The distance distribution profiles between Ca atoms of residues L44 and V55 are shown in (a). The effect of T54I mutation in 
bridging the two adjacent sheets is shown in (b,c). 



at position T54 has been established (Welsch et al., 
2008). In the structure of HCV-lb NS3 protease, T54 
is located at the very end of a (3-strand (Figure 4(c)), 
which belongs to an antiparallel P-sheet (Welsch et al., 
2008). The hydroxyl group of the T54 side chain is 
involved in the formation of two H-bonds with the 
adjacent residues V55 and L44 in the neighboring 
antiparallel (3-strand (Figure 4(c)). The distance 
between the back-bone H-bond donor and acceptor, 
L44 and V55, at 4.60 A is too large to be bridged by 
a single H-bond (Welsch et al, 2008). The two H- 
bonds from the T54 residue side chain bridge the two 
strands and thereby stabilize the local P-sheet confor- 
mation (Figure 4(c)). 

The sequence of HCV-4a is different from HCV-lb 
at residue 54 (Figure 1(a) and (d)). The absence of a 
hydroxyl group in T541 is expected to have an impact 
on the local P-sheet conformation in HCV-4a protease 
(Figure 4(b)). To assess the dynamics consequences of 
that mutation, we used the distance between the Ca 
atoms of amino acids L44 and V55 as an indicator of 
the conformational changes. We observed a broader 
distance distribution (~50% increase) and an upshift in 
peak value(-0.3A) between HCV-4a and HCV-lb 
(Figure 4(a)). Thus, the strands bearing L44 and V55 
are expected to spend more time closer to each other 
in HCV-lb than with HCV-4a. 

As H57 is adjacent to V55, any loss of stability in 
V55, due to the lack of H-bond with 154, may directly 
impact the dynamics of the catalytic histidine (Figure 4 
(b)). Thus, the larger oscillation of the V55 bearing 
strand in HCV-4a when coupled with the H57's large 
RMSD, points to a histidine that is less capable to 
perform its catalytic function. This could be an impor- 
tant contribution not only to the higher drug resistivity, 



but also to the decreased protease activity of HCV-4a 
compared with HCV-lb as mentioned in Massariol, 
Zhao, Marquis, Thibeault, and White (2010). 

Thirty variability sites exist between the sequences of 
HCV-lb and HCV-4a (Figure 1(d)) and it is unlikely that 
T541, alone, is responsible for the differences in dynam- 
ics at the catalytic site. In other words, the genotypic 
variations express their effects in combination and not as 
individual mutations. Therefore, our data is relevant only 
when all these genotypic variations are present in the 
structure. 

Taken together, the data presented here suggest that 
the structural dynamics within the catalytic triad are sig- 
nificantly different between HCV-lb and HCV-4a NS3 
proteases. The effect of sequence variability in the HCV- 
4a protease seems to extend beyond the rigid structure to 
render the molecular dynamics less accommodating to 
natural substrates and, possibly, potential drugs. It is evi- 
dent that there are numerous details underpinning the 
causes of the genotype-based drug resistance. However, 
the changes in the dynamic behavior at the catalytic 
triad, between drug responsive and drug resistant HCV 
genotypes, provide a new avenue of inquiry with sugges- 
tive results. The final confirmation awaits the crystal 
structures of HCV-4a protease, with and without bound 
drugs, along with computationally intensive analyses of 
the protease-drug complexes, and the interplay between 
sequence variability, dynamic behavior, and drug 
binding. 

Supplementary material 

The supplementary material for this paper is available 
online at http://dx.doi.10.1080/07391102.2013.800001. 
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